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Abstract. Implicit and explicit raytracing-photoionisation algorithms have been im- 
plemented in the author's radiation-magnetohydrodynamics code. The algorithms are 
described briefly and their efficiency and parallel scaling are investigated. The implicit 
algorithm is more efficient for calculations where ionisation fronts have very supersonic 
velocities, and the explicit algorithm is favoured in the opposite limit because of its bet- 
ter parallel scaling. The implicit method is used to investigate the effects of initially 
uniform magnetic fields on the formation and evolution of dense pillars and cometary 
globules at the boundaries of HII regions. It is shown that for weak and medium field 
strengths an initially perpendicular field is swept into alignment with the pillar during its 
dynamical evolution, matching magnetic field observations of the 'Pillars of Creation' 
in Ml 6. A strong perpendicular magnetic field remains in its initial configuration and 
also confines the photoevaporation flow into a bar-shaped, dense, ionised ribbon which 
partially shields the ionisation front. 



1. Introduction 



Motivated by observations o f ongoing star form ation within pillars and globules at the 
borders of H n regions (e.g. iHester et all 1 19961) . and by increasing evidence that star 
formation propaga tes outwards from young clusters sequentially dLee & Chenl 120071 : 
Smith et al. I l2010bh . models for the formation and evolution of theses structures have 
been investigated in increasing detail recently. These are highly non-linear structures 
forming in the presence of intense radiation fields, so multidimensional numerical sim- 
ulations including hydrodynamics (HD) with photoionisation (R-HD) are required, ide- 
ally also including self-gravity and magnetohydrody namics (MH D). Low resolution 2D 
calculations were already performed in the 1980s (lYorkd[l986h : ionisation front insta- 
bilities were shown to produce structures resembling pillars during the expansion of 
H ii regions from the ultra-compac t phase dGarcfa-Segura & Francdll996h . also shown 
more recently in 3D simulations dWhalen & NormarJl2008h . Shadowing due to pre- 
existing o verdensities was als o shown to generate pillar-like structures in R-HD sim- 
ulations ( Wil liams et alJl200l | ): the effectiveness of thi s process in a turbulent de nsity 
field was demonstrated by iMellema et al J ([2006al) and iGrits chneder et all (1201 Oh. and 
in the presence of magnetic fields by I Arthur et al.1 (1201 lh . In lMackey & Liml (1201 Oh we 
studied the formation of massive pillars due to shadowing of randomly placed dense 
clumps with 3D R-HD simulations, showing that this simple model could reproduce 
the main features of optical images, the morphology of molecular emission, and the 
line-of-sight v elocity profiles of the M 16 pillars obtained from molecular line emission 
(|Poundl [l998). In this contribution simulations performed by the author and Dr. An- 
drew J. Lim are described which include MHD with non-equilibrium photoionisation 
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and raytracing (R-MHD), with the aim of determining t he ef fects of magnetic fields on 
the pillar formati on process stud i ed in Mackey & Liml (l2010h . This work is described 
in more detail in lMackey & Liml (|2011) and is c o mplem entary to R-MHD simulations 
of global H ii region evolution by I Arthur et all (1201 lh and photoionisation of single 
clumps by Henney et al. (2009). In the following section the R-MHD code used for the 
calculations is described and its scaling behaviour on parallel supercomputers is dis- 
cussed for implicit and explicit raytracing-photoionisation algorithms. In section [3] the 
R-MHD simulations are described and the results obtained are summarised. 



2. Code Description and Parallel Scaling 



In our R-MHD code dMackey & Limll2010l . 1201 lh the equations of ideal compressible 



MHD are solved on a uniform grid using a directionally-unsplit, second order (in time 
and space), finite volume formulation with a Roe-type Riemann solver in conserved 
variables. A tracer variable is used for the ion fraction of H; microphysical (radia- 
tive and collisional) heating and cooling processes also provide a source term to the 
energy equation and are included by operator splitting. The column density from a ra- 
diation source to a given grid cell is calculated using a short characteristics raytracing 
module; diffuse radiation is treated approximately by the on-the-spot approximation. 
Photoionisation heating and recombination cooling are calculated explicitly for H us- 
ing its non-equili brium ion fraction; coo ling due to other elements is approximate and 
uses model C2 in lMackey & Liml (l2010h . The code is parallelised using MPI and was 



run on 128 cores for the calculations in section [3j 

An implicit raytracing algorithm w as used for the simula tions presented here, 
based closely on the C 2 -ray algorithm of Mell ema et al.1 (l2006bh but with some mod- 



ifications (referred to here as Algorithm 1 or just Algl). This updates microphysical 
quantities (internal energy, ion fraction) as rays are being traced outwards from the 
source, using time-averaged attenuation fractions through cells and thus allowing ioni- 
sation fronts to cross many optically thick grid cells in a single raytracing step without 
loss of photon conservation. The only serious disadvantage of this integration scheme 
is that it is difficult to parallelise efficiently on non-spherical grids: raytracing outwards 
from a source must proceed sequentially from one processor's subdomain to its neigh- 
bour and so on. Since the microphysical quantities are integrated during raytracing this 
places a significant fraction of the total computation in a poorly- scaling algorithm. In 
the limit of many subdomains the parallel scaling of this algorithm is simple to deduce: 
the radial direction must be calculated (at least partially) in serial and only the perpen- 
dicular direction(s) can be efficiently parallelised. In 2D, therefore, the strong scaling 
of the raytracing algorithm on N processes will have runtime t s { m oc Af" 1/2 , and in 3D 
the scaling will be t s [ m oc N~ 213 . This scaling holds for the short characteristics tracer 
(and probably also for other ray-splitting algorithms). 

An explic it raytracing algorithm ha s recently been added to the code, similar to 
the method of IWhalen & Normianl (l2006h . hereafter denoted Alg3. Here only the col- 
umn densities are calculated in the raytracing step and the microphysics integration is 
performed fully in parallel. A second order integration, using the same approach as 
for the hydrodynamics, proved much more efficient than the usual first order method. 
A raytracing is performed at the start of the timestep; everything is integrated forwards 
half a timestep; another raytracing is performed using this time-centred density (and ion 
fraction) field; these column densities are then used to integrate the full timestep from 
the initial to time-advanced values. With this algorithm only 8 raytracings (4 ti mesteps) 
are required to photoionise an optically thick cell with an error of less than 1% (iMackeyl 
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Figure 1. Scaling comparison between Algorithms 1 and 3 on the distributed 
memory cluster JUROPA. Left: Scaling of a 2D static problem as discussed in the 
text. Right: Scaling for a 2D calculation with photoionisation and a 2000 km s" 1 
stellar wind, as discussed in the text. Timing was measured with a microsecond 
timer in the main simulation loop, not including time spent in initialisation or data 
output. 



20111 ): the microphysics timestep limit used is At = 0.25/i, where x is the time rate of 



change of the H + fraction, x. 

The total runtime for the code in 2D simulations (cylindrical coordinates [R, z]) as 
a function of number of processors N is shown in Fig. Q] for two different calculations. 
The first simulation is of the expansion of an H n region into a uniform static medium 
covering the range [R,z] e [0, 2] pc with a constant gas density of nu = 100 cm" 3 and 
initial temperature T = 50K. The star at the origin has monochromatic ionising photon 
luminosity N = 10 48 s" 1 . The grid contains 256 2 cells and the simulation was run for 
10 recombination times (f rec , here 3.861 x 10 10 s) on 8 - 128 processors (Fig. [TJ left 
panel). The better scaling properties of the explicit Alg3 are evident, but the overall 
efficiency of Algl is so much greater that it remains the faster algorithm. 

In the second simulation the star at the origin has luminosity N = 3 x 10 48 s" 1 and 
a fast stellar wind is injected in a 10 cell radius around the origin with mass-loss rate 
M = 10" 7 M yr _1 and velocity v w = 2000 km s" 1 . The model is run with 512 x 256 
grid cells and physical domain z e [-1.28, 1.28] x 10 18 cm, R e [0, 1.28] x 10 18 cm, 
with ambient density nu = 3000 cm" 3 . Here the overall timestep is set by the Courant 
condition on the 2000 km s" 1 wind region and not by the microphysics. The simulation 
is run for 10 kyr, or ^ 245^ rec . The right panel of Fig. Q] shows that Alg3 is more 
efficient than Algl for this problem on 32-128 cores. This new algorithm will be used 
for calculations of the circumstellar medium around evolving massive stars. 



3. R-MHD Simulation Results 

We have performed a number of 3D R-MHD simulat ions of the formation a nd evolu- 
tion of pillars and globules, presented in more detail in lMackey & Liml (1201 lh . Here we 



present results from simulations involving the photoionisation of three dense clumps of 
gas embedded in a uniform ambient medium and exposed to ionising radiation from a 
point source. The ambient gas density is tin = 200 cm" 3 , and the clumps have peak 
densities 500x greater with Gaussian density profiles [p oc exp(-r 2 /2rQ)] of width 
ro = 0.09 pc. The 3D simulation domain is x e [1.5,6.0] pc, {y,z} e [-1.5, 1.5] pc 
resolved by 384 x 256 2 grid cells; the radiation source is off the grid at the origin with 
monochromatic photon luminosity N = 2 x 10 50 s _1 . The three clumps are almost 
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Figure 2. Projected gas density (x-z plane) for models R2, R5 and R8 on the 
logarithmic column density scale indicated. The radiation source is located above 
the simulation domain at [0,0,0]. Projected magnetic field orientation is overlaid 
(weighted by gas density in the line-of-sight integral, and recovered from Stokes Q 
and U parameters). The left panel shows R2 at t = 25 kyr, the other three show R2, 
R5, and R8 (from left to right) at t = 250 kyr. 

collinear with positions (in pc) [2.3, 0, 0], [2.75, 0, 0.12] and [3.2, 0, -0.12]. Three sim- 
ulations are presented here: R2 has an initial magnetic field B = [1.8, 1.8, 17.7] yt/G, R5 
has B = [0,0, 53.2] yuG, and R8 has B = [14.2,7.1, 159]yuG. All three models have a 
field largely perpendicular to the radiation propagation direction. The weak B field in 
R2 means that the re sults closely follow purely hydrodynamic evolution (see fig. 7 in 
Mackey & Lim|[20Tol . describing model 17), whereas R8 is largely magnetically dom- 
inated; R5 lies between the two extremes. All models were evolved to at least 400 
kyr. 

The projected neutral hydrogen column density is shown in Fig. [2| at times / = 
25 kyr (for R2) and 250 kyr (R2, R5, R8), with projected magnetic field orientation 
indicated by the overlaid black lines. At 25 kyr all three models are effectively identical 
- an R-type ionisation front has ionised most of the low density gas on the grid except 
for the shadowed region and the three clumps are beginning to be compressed by shocks 
generated by the high pressure ionised gas. By 250 kyr, however, clear differences have 
emerged: the weak field model R2 behaves almost identically to a purely hydrodynamic 
simulation but the strong field model R8 is very diffeent. In dense neutral gas in R2 the 
initially perpendicular magnetic field has been swept into alignment with the long axis 
of the dense pillar. The projected B field in the ionised gas around the head of the 
pillar follows the radial photoevaporation flow generated at the ionisation front; this is 
also seen in model R5 but is absent in R8. The alignment of the magnetic field with 
the long axis of the pillar is also weaker in R5 and the field orientation is essentially 
unchanged in R8. The projected density structure is rather different in R8 also, in that 
the clump furthest from the source has stayed largely intact and separate from the two 
closer to the source (which have merged). In R2 an R5 the furthest clump is no longer 
distinguishable from the general structure of the pillar. 

R8 is also different in that it has significant neutral low-column-density gas (dark 
blue). This can be understood by looking at the ionised gas emission in Fig. [3] which 
shows that in R8 there is a lot of dense ionised gas ahead of the pillar. This provides 
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Figure 3. Simulations R2 (above), R5 (centre) and R8 (below) at t = 250 kyr, 
shown here in Ha emission (arbitrary logarithmic intensity scale). The left sequence 
shows projections along the z axis (parallel to initial magnetic field orientation), 
while the right hand side shows projections along the y axis, perpendicular to the 
field and the same projection as in Fig. [2J 



significant optical depth and allows gas to recombine in its shadow around the pillar. 
This gas has been organised by the strong magnetic field into a bar-shaped or filamen- 
tary structure. The photoevaporation flow is clearly seen in these figures for models 
R2 and R5, but for R8 it is confined to a much smaller volume, set by the radius at 
which the ram pressure of the flow equals the magnetic pressure in the ionised gas, at 
which point an isothermal reverse shock forms. The shocked ionised gas then disperses 
subsonically, but to first order it is only free to move along field lines so it stretches into 
the filamentary structure seen in this figure. 



4. Conclusions 

Tests of implicit and explicit raytracing algorithms show that for simulations without 
rapidly moving (v » lOkms" 1 ) or high temperature (T » 10 4 K) gas the implicit 
method is always more efficient on the uniform grid used here because of the large 
timesteps it allows. By contrast, for a model with a fast stellar wind (v = 2000 km s" 1 ) 
the superior parallel scaling of the explicit algorithm makes it more efficient. 

Our R-MHD photoionisation simulations show that both radiation-driven implo- 
sion and acceleration of clumps by the rocket effect tend to align the mag netic field with 
the ra diation propagation direction in dense neutral gas (suggested by Sugitani et al. 



20071 for the case of M 16). The effectiveness of this alignment is dependent on the 
initial field strength: with similar gas densities and pressures to those in M 16, these 
simulations show that a magnetic field of 160yt/G is sufficient to prevent any significant 
field reorientation. By interpolating between the simulation results, an ambient field of 
| B| < 50 A/G is requi red to explain the observed field configuration in the M 16 pillars 



(Sugi tani et al.ll2007b if the pillars formed via the mechanism we are modelling. If the 



field were stronger than this, additional signatures would also be seen in the ionised gas 
emission tracing the large-scale field orientation. Such features appear to be present in 
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Ha images of the pillar near N GC 6357 dBohigas et alJl2004l) and in a pillar in the Ca 



rina nebula ( Smit h et aDl2010al . Fig. 1). These results show that magnetic field strengths 



in star-forming regions can in principle be significantly constrained by the morphology 
of structures which form at the borders of HII regions. 
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